############################################
## author:    Robert A. Huber
## contact:   robert.huber@sbg.ac.at
## file name: beliefs_regression1.R
## Context:   Environmental Politics Paper on Beliefs
## started:   2018-10-05
## Summary:   Robustness Checks for analyses of the first conjoint
##############################################
 
df_conj_out1$eff <- abs(df_conj_out1$eff - 8)
df_conj_out1$edu <- factor(df_conj_out1$edu)

# Explanations of Beliefs ####

# Regresison Analyses ####

m1c <- lmer(choice ~ eff + fai + int + (1 | id), data = df_conj_out1)

m1cCL <- lm.cluster(formula = choice ~ eff + fai + int, data = df_conj_out1, cluster = "id")

m1cFE <- lm(choice ~ eff + fai + int + factor(id), data = df_conj_out1)

m2c <- lmer(choice ~ eff * fai * int + (1 | id), data = df_conj_out1)

m2cCL <- lm.cluster(formula = choice ~ eff * fai * int, data = df_conj_out1, cluster = "id")

m2cFE <- lm(choice ~ eff * fai * int + factor(id), data = df_conj_out1)

m1r <- lmer(rate ~ eff + fai + int + frame + age + gender + edu + mp + uar + leftright + (1 | id), data = df_conj_out1)

m1rCL <- lm.cluster(formula = rate ~ eff + fai + int + frame + age + gender + edu + mp + uar + leftright, data = df_conj_out1, cluster = "id")

m1rFE <- lm(rate ~ eff + fai + int + factor(id), data = df_conj_out1)

m2r <- lmer(rate ~ eff*fai*int + frame + age + gender + edu + mp + uar + leftright + (1 | id), data = df_conj_out1)

m2rCL <- lm.cluster(formula = rate ~ eff * fai * int + frame + age + gender + edu + mp + uar + leftright, data = df_conj_out1, cluster = "id")

m2rFE <- lm(rate ~ eff * fai * int + factor(id), data = df_conj_out1)

m3r <- lmer(rate ~ eff*frame + fai*frame + int*frame + frame + age + gender + edu + mp + uar + leftright + (1 | id), data = df_conj_out1)

m3rCL <- lm.cluster(formula = rate ~ eff*frame + fai*frame + int*frame + frame + age + gender + edu + mp + uar + leftright, data = df_conj_out1, cluster = "id")

m3rFE <- lm(rate ~ eff*frame + fai*frame + int*frame + factor(id), data = df_conj_out1)

texreg(list(m1c, m2c, m1r, m2r, m3r),
       #longtable = T,
       sideways = T,
       leading.zero = T,
       digits = 2,
       single.row = T,
       reorder.coef = c(2:4,9:15,16:17, 5:8,18:21, 1),
       custom.model.names = c("Choice: Baseline", "Choice: 3way Interaction",
                              "Rating: Baseline", "Rating: 3way Interaction", 
                              "Rating by Frame"),
       custom.coef.names = c("Intercept", "Effectiveness", "Fairness",
                             "Intrusiveness", "Effectiveness x Fairness",
                             "Effectiveness x Intrusiveness", "Intrusiveness x Fairness",
                             "Effectiveness x Intrusiveness x Fairness",
                             "Frame (Ref = EV)", "Age",
                             "Gender (Ref = Female)", "Education Other",
                             "Education - Secondary", "Education - Tertiary",
                             "Mobility Profile", "Living Situation - Agglomeration",
                             "Living Situation - Rural", "Political Ideology", "Effectiveness x Frame",
                             "Fairness x Frame", "Intrusiveness x Frame"),
       custom.note = "%stars. Entries are unstandardised regression coefficients from a linear mixed effects regression.",
       custom.gof.names = c(NA, NA, NA, "Number of observations", "Number of individuals", NA, NA))

texreg(list(m1cFE, m2cFE, m1rFE, m2rFE, m3rFE),
       leading.zero = T,
       digits = 2,
       caption = "Perceived Policy Consequences and Policy Support using Fixed Effects",
       caption.above = T,
       sideways = T,
       single.row = T,
       omit.coef = "id",
       reorder.coef = c(2:12,1),
       custom.model.names = c("Choice: Baseline", "Choice: 3way Interaction",
                              "Rating: Baseline", "Rating: 3way Interaction", 
                              "Rating by Frame"),
       custom.coef.names = c("Intercept", "Effectiveness", "Fairness",
                             "Intrusiveness", "Effectiveness x Fairness",
                             "Effectiveness x Intrusiveness", "Intrusiveness x Fairness",
                             "Effectiveness x Intrusiveness x Fairness",
                             "Frame (Ref = EV)", "Effectiveness x Frame",
                             "Fairness x Frame", "Intrusiveness x Frame"),
       custom.note = "%stars. Entries are unstandardised regression coefficients from a linear regression with individual fixed effects.",
       custom.gof.names = c(NA, NA, "Number of observations", "RMSE"))

texreg(list(m1c, m2c, m1r, m2r, m3r),
       leading.zero = T,
       digits = 2,
       single.row = T,
       caption = "Perceived Policy Consequences and Policy Support using Clustered Standard Errors",
       caption.above = T,
       sideways = T,
       override.coef = list(summary(m1cCL)[,1], summary(m2cCL)[,1], summary(m1rCL)[,1], summary(m2rCL)[,1], summary(m3rCL)[,1]),
       override.se = list(summary(m1cCL)[,2], summary(m2cCL)[,2], summary(m1rCL)[,2], summary(m2rCL)[,2], summary(m3rCL)[,2]),
       override.pvalues = list(summary(m1cCL)[,4], summary(m2cCL)[,4], summary(m1rCL)[,4], summary(m2rCL)[,4], summary(m3rCL)[,4]),
       reorder.coef = c(2:4,9:15,16:17, 5:8,18:21, 1),
       custom.model.names = c("Choice: Baseline", "Choice: 3way Interaction",
                              "Rating: Baseline", "Rating: 3way Interaction", 
                              "Rating by Frame"),
       custom.coef.names = c("Intercept", "Effectiveness", "Fairness",
                             "Intrusiveness", "Effectiveness x Fairness",
                             "Effectiveness x Intrusiveness", "Intrusiveness x Fairness",
                             "Effectiveness x Intrusiveness x Fairness",
                             "Frame (Ref = EV)", "Age",
                             "Gender (Ref = Female)", "Education Other",
                             "Education - Secondary", "Education - Tertiary",
                             "Mobility Profile", "Living Situation - Agglomeration",
                             "Living Situation - Rural", "Political Ideology", "Effectiveness x Frame",
                             "Fairness x Frame", "Intrusiveness x Frame"),
       custom.note = "%stars. Entries are unstandardised regression coefficients from a linear regression with clustered standard errors.",
       custom.gof.names = c(NA, NA, NA, "Number of observations", "Number of individuals", NA, NA))


# Test 3w Interaction using Effects Package -------------------------------

#rating Effectiveness
fit.effR <- effect("eff:fai:int", m2r, xlevels = list(eff=c(1:7), 
                                                      fai=c(mean(df_conj_out1$fai)-sd(df_conj_out1$fai), mean(df_conj_out1$fai)+sd(df_conj_out1$fai)),
                                                      int=c(mean(df_conj_out1$int)-sd(df_conj_out1$int), mean(df_conj_out1$int)+sd(df_conj_out1$int))),
                   x.var="eff", confidence.level = 0.95)

ratingE <- data.frame(fit.effR)

ratingE$Y <- ifelse(ratingE$int == min(ratingE$int), "Low Y", "High Y")
ratingE$Z <- ifelse(ratingE$fai == min(ratingE$fai), "Low Z", "High Z")
ratingE$YZ <- factor(paste(ratingE$Y, ratingE$Z, sep=" / "))

ratingE$Y_lab <- ifelse(ratingE$int == min(ratingE$int), "Low Intrusiveness", "High Intrusiveness")
ratingE$Z_lab <- ifelse(ratingE$fai == min(ratingE$fai), "Low Fairness", "High Fairness")
ratingE$YZ_lab <- factor(paste(ratingE$Y_lab, ratingE$Z_lab, sep=" / "))

#ratingIntrusiveness
fit.effR <- effect("eff:fai:int", m2r, xlevels = list(int=c(1:7), 
                                                      fai=c(mean(df_conj_out1$fai)-sd(df_conj_out1$fai), mean(df_conj_out1$fai)+sd(df_conj_out1$fai)),
                                                      eff=c(mean(df_conj_out1$eff)-sd(df_conj_out1$eff), mean(df_conj_out1$eff)+sd(df_conj_out1$eff))),
                   x.var="int", confidence.level = 0.95)

ratingI <- data.frame(fit.effR)

ratingI$Y <- ifelse(ratingI$eff == min(ratingI$eff), "Low Y", "High Y")
ratingI$Z <- ifelse(ratingI$fai == min(ratingI$fai), "Low Z", "High Z")
ratingI$YZ <- factor(paste(ratingI$Y, ratingI$Z, sep=" / "))

ratingI$Y_lab <- ifelse(ratingI$eff == min(ratingI$eff), "Low Effectiveness", "High Effectiveness")
ratingI$Z_lab <- ifelse(ratingI$fai == min(ratingI$fai), "Low Fairness", "High Fairness")
ratingI$YZ_lab <- factor(paste(ratingI$Y_lab, ratingI$Z_lab, sep=" / "))

#ratingIntrusiveness
fit.effR <- effect("eff:fai:int", m2r, xlevels = list(fai=c(1:7), 
                                                      int=c(mean(df_conj_out1$int)-sd(df_conj_out1$int), mean(df_conj_out1$int)+sd(df_conj_out1$int)),
                                                      eff=c(mean(df_conj_out1$eff)-sd(df_conj_out1$eff), mean(df_conj_out1$eff)+sd(df_conj_out1$eff))),
                   x.var="fai", confidence.level = 0.95)

ratingF <- data.frame(fit.effR)

ratingF$Y <- ifelse(ratingF$eff == min(ratingF$eff), "Low Y", "High Y")
ratingF$Z <- ifelse(ratingF$int == min(ratingF$int), "Low Z", "High Z")
ratingF$YZ <- factor(paste(ratingF$Y, ratingF$Z, sep=" / "))

ratingF$Y_lab <- ifelse(ratingF$eff == min(ratingF$eff), "Low Effectiveness", "High Effectiveness")
ratingF$Z_lab <- ifelse(ratingF$int == min(ratingF$int), "Low Intrusiveness", "High Intrusiveness")
ratingF$YZ_lab <- factor(paste(ratingF$Y_lab, ratingF$Z_lab, sep=" / "))

#Choice Effectiveness
fit.effR <- effect("eff:fai:int", m2c, xlevels = list(eff=c(1:7), 
                                                      fai=c(mean(df_conj_out1$fai)-sd(df_conj_out1$fai), mean(df_conj_out1$fai)+sd(df_conj_out1$fai)),
                                                      int=c(mean(df_conj_out1$int)-sd(df_conj_out1$int), mean(df_conj_out1$int)+sd(df_conj_out1$int))),
                   x.var="eff", confidence.level = 0.95)

choiceE <- data.frame(fit.effR)

choiceE$Y <- ifelse(choiceE$int == min(choiceE$int), "Low Y", "High Y")
choiceE$Z <- ifelse(choiceE$fai == min(choiceE$fai), "Low Z", "High Z")
choiceE$YZ <- factor(paste(choiceE$Y, choiceE$Z, sep=" / "))

choiceE$Y_lab <- ifelse(choiceE$int == min(choiceE$int), "Low Intrusiveness", "High Intrusiveness")
choiceE$Z_lab <- ifelse(choiceE$fai == min(choiceE$fai), "Low Fairness", "High Fairness")
choiceE$YZ_lab <- factor(paste(choiceE$Y_lab, choiceE$Z_lab, sep=" / "))

#Choice Intrusiveness
fit.effR <- effect("eff:fai:int", m2c, xlevels = list(int=c(1:7), 
                                                      fai=c(mean(df_conj_out1$fai)-sd(df_conj_out1$fai), mean(df_conj_out1$fai)+sd(df_conj_out1$fai)),
                                                      eff=c(mean(df_conj_out1$eff)-sd(df_conj_out1$eff), mean(df_conj_out1$eff)+sd(df_conj_out1$eff))),
                   x.var="int", confidence.level = 0.95)

choiceI <- data.frame(fit.effR)

choiceI$Y <- ifelse(choiceI$eff == min(choiceI$eff), "Low Y", "High Y")
choiceI$Z <- ifelse(choiceI$fai == min(choiceI$fai), "Low Z", "High Z")
choiceI$YZ <- factor(paste(choiceI$Y, choiceI$Z, sep=" / "))

choiceI$Y_lab <- ifelse(choiceI$eff == min(choiceI$eff), "Low Effectiveness", "High Effectiveness")
choiceI$Z_lab <- ifelse(choiceI$fai == min(choiceI$fai), "Low Fairness", "High Fairness")
choiceI$YZ_lab <- factor(paste(choiceI$Y_lab, choiceI$Z_lab, sep=" / "))

#Choice Intrusiveness
fit.effR <- effect("eff:fai:int", m2c, xlevels = list(fai=c(1:7), 
                                                      int=c(mean(df_conj_out1$int)-sd(df_conj_out1$int), mean(df_conj_out1$int)+sd(df_conj_out1$int)),
                                                      eff=c(mean(df_conj_out1$eff)-sd(df_conj_out1$eff), mean(df_conj_out1$eff)+sd(df_conj_out1$eff))),
                   x.var="fai", confidence.level = 0.95)

choiceF <- data.frame(fit.effR)

choiceF$Y <- ifelse(choiceF$eff == min(choiceF$eff), "Low Y", "High Y")
choiceF$Z <- ifelse(choiceF$int == min(choiceF$int), "Low Z", "High Z")
choiceF$YZ <- factor(paste(choiceF$Y, choiceF$Z, sep=" / "))

choiceF$Y_lab <- ifelse(choiceF$eff == min(choiceF$eff), "Low Effectiveness", "High Effectiveness")
choiceF$Z_lab <- ifelse(choiceF$int == min(choiceF$int), "Low Intrusiveness", "High Intrusiveness")
choiceF$YZ_lab <- factor(paste(choiceF$Y_lab, choiceF$Z_lab, sep=" / "))

df_3w <- rbind(ratingE, ratingI, ratingF, choiceE, choiceI, choiceF)
df_3w$task <- rep(c("Rating", "Choice"), each = nrow(df_3w)/2)
df_3w$task <- factor(df_3w$task, levels = c("Choice", "Rating"))
df_3w$IV <- rep(rep(c("Effectiveness", "Intrusiveness", "Fairness"), each = nrow(df_3w)/6),2)
df_3w$IV <- factor(df_3w$IV, levels = c("Effectiveness", "Intrusiveness", "Fairness"))
df_3w$X <- ifelse(df_3w$IV == "Effectiveness", df_3w$eff,
                  ifelse(df_3w$IV == "Intrusiveness", df_3w$int,
                         ifelse(df_3w$IV == "Fairness", df_3w$fai, NA)))

#What is Y and Z -> text_annotation

figureA5 <- ggplot(subset(df_3w, df_3w$IV == "Effectiveness"), aes(x=X, y=fit)) +
  geom_line(aes(colour = YZ_lab), lwd = 1.05) +
  geom_line(aes(y=lower,colour = YZ_lab), lty= "dashed", lwd = .5) +
  geom_line(aes(y=upper,colour = YZ_lab), lty= "dashed", lwd = .5) +
  geom_point(aes(shape = YZ_lab, colour = YZ_lab), size = 2.5) +
  theme_bw() +
  facet_grid(task~.) +
  ylab("Prediction") + 
  xlab("Belief") +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_color_manual(name = "Legend", values = rep(c(1:4), 3))+
  scale_shape_manual(name = "Legend", values = rep(c(15:19), 3)) +
  theme(text = element_text(size=20), axis.text.x = (element_text(angle=25, hjust=1)))
figureA5

figureA6 <- ggplot(subset(df_3w, df_3w$IV == "Intrusiveness"), aes(x=X, y=fit)) +
  geom_line(aes(colour = YZ_lab), lwd = 1.05) +
  geom_line(aes(y=lower,colour = YZ_lab), lty= "dashed", lwd = .5) +
  geom_line(aes(y=upper,colour = YZ_lab), lty= "dashed", lwd = .5) +
  geom_point(aes(shape = YZ_lab, colour = YZ_lab), size = 2.5) +
  theme_bw() +
  facet_grid(task~.) +
  ylab("Prediction") + 
  xlab("Belief") +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_color_manual(name = "Legend", values = rep(c(1:4), 3))+
  scale_shape_manual(name = "Legend", values = rep(c(15:19), 3))+
  theme(text = element_text(size=20), axis.text.x = (element_text(angle=25, hjust=1)))
figureA6

figureA7 <- ggplot(subset(df_3w, df_3w$IV == "Fairness"), aes(x=X, y=fit)) +
  geom_line(aes(colour = YZ_lab), lwd = 1.05) +
  geom_line(aes(y=lower,colour = YZ_lab), lty= "dashed", lwd = .5) +
  geom_line(aes(y=upper,colour = YZ_lab), lty= "dashed", lwd = .5) +
  geom_point(aes(shape = YZ_lab, colour = YZ_lab), size = 2.5) +
  theme_bw() +
  facet_grid(task~.) +
  ylab("Prediction") + 
  xlab("Belief") +
  theme(strip.background = element_rect(fill = 'white')) +
  scale_color_manual(name = "Legend", values = rep(c(1:4), 3))+
  scale_shape_manual(name = "Legend", values = rep(c(15:19), 3))+
  theme(text = element_text(size=20), axis.text.x = (element_text(angle=25, hjust=1)))
figureA7

rm(list=setdiff(ls(), c("ri", "df_conj_out1", "df_conj_out2")))
